
#include "gausssolver.h"
#include <cmath>

#include <omp.h>
#include <iostream>

void GaussSolver::solve() {
	if (matrix.getDimension() != freeTerms.size()) {
		setState(ssNotSolvable);
		return;
	}

	int sz = freeTerms.size();
	Xs.resize(sz);
	for (int i = 1; i <= sz; i++) {
		if (matrix(i, i) == 0) {
			int l = i;
			while (l <= sz && matrix(l, i) == 0) l++;
			if (l <= sz) {
				#pragma omp parallel
				{
					//std::cout << omp_get_num_threads() << std::endl;
					#pragma omp for
					for (int k = 1; k <= sz; k++) {
						#pragma omp critical
						{
							double t = matrix(i, k);
							matrix(i, k) = matrix(l, k);
							matrix(l, k) = t;
						}
					}
					#pragma omp critical
					{
						double t = freeTerms[i-1];
						freeTerms[i-1] = freeTerms[l-1];
						freeTerms[l-1] = t;
					}
				}
			}
		}
		if (matrix(i, i) != 0) {
			for (int j = i+1; j <= sz; j++) {
				#pragma omp parallel
				{
					//std::cout << omp_get_num_threads() << std::endl;
					#pragma omp for
					for (int k = i+1; k <= sz; k++) {
						#pragma omp critical
						matrix(j, k) -= matrix(j, i) / matrix(i, i) * matrix(i, k);
					}
					#pragma omp critical
					{
						freeTerms[j-1] -= matrix(j, i) / matrix(i, i) * freeTerms[i-1];
						matrix(j, i) = 0;
					}
				}
			}
		} else {
			setState(ssNotSolvable);
			return;
		}
	}
	for (int i = sz; i > 0; i--) {
		#pragma omp parallel
		{
			//std::cout << omp_get_num_threads() << std::endl;
			#pragma omp critical
			Xs[i-1] = freeTerms[i-1];
			#pragma omp for
			for (int k = i+1; k <= sz; k++) {
				#pragma omp critical
				Xs[i-1] -= matrix(i, k) * Xs[k-1];
			}
			#pragma omp critical
			Xs[i-1] /= matrix(i, i);
		}
	}
}
